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Abstract. We report the results of a detailed study of the spectral properties of Laplace and 
Stokes operators, modified with a volume penalization term designed to approximate Dirichlet condi- 
tions in the limit when a penalization parameter, r], tends to zero. The eigenvalues and eigenfunctions 
are determined either analytically or numerically as functions of r], both in the continuous case and af- 
ter applying Fourier or finite difference discretization schemes. For fixed r\, we find that only the part 
of the spectrum corresponding to eigenvalues A < rf~ 1 approaches Dirichlet boundary conditions, 
while the remainder of the spectrum is made of uncontrolled, spurious wall modes. The penalization 
error for the controlled eigenfunctions is estimated as a function of 77 and A. Surprisingly, in the 
Stokes case, we show that the eigenfunctions approximately satisfy, with a precision 0{rf), Navier 
slip boundary conditions with slip length equal to ^/fj. Moreover, for a given discretization, we show 
that there exists a value of 7], corresponding to a balance between penalization and discretization 
errors, below which no further gain in precision is achieved. These results shed light on the behavior 
of volume penalization schemes when solving the Navier-Stokes equations, outline the limitations of 
the method, and give indications on how to choose the penalization parameter in practical cases. 
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1. Introduction. Penalization approaches, as reviewed for example in Ref. [221 
118] . are nowadays commonly employed to solve boundary or initial-boundary value 
problems, as encountered for example in computational fluid dynamics. They consist 
in embedding the original, possibly complex spatial domain inside a bigger domain 
having a simpler geometry, for example a torus, while keeping the boundary conditions 
approximately enforced thanks to new terms that are added to the equations. One 
particular example is the volume penalization method which, inspired by the physical 
intuition that a solid wall in contact with a fluid is similar to a vanishingly porous 
medium, uses the Brinkman-Darcy drag force [fH [2] as penalization term. Due to 
this analogy, the method is sometimes called Brinkman penalization. It has been 
mathematically justified in academic cases [HEj, an d its use has rapidly spread to 
various domains of scientific computing. 

The main advantage of such penalized equations is that they can be discretized 
independently of the geometry of the original problem, as the latter has been encoded 
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into the penalization terms. Such a simplification permits a massive reduction in solver 
development time, since it avoids the issues associated to the design and management 
of the grid, allowing for example the use of simple spectral solvers in Cartesian ge- 
ometries [TTJ [TU [21] . The gain becomes even more substantial when the geometry is 
time-dependent, as in the case of moving obstacles [TBI [15], or when fluid-structure 
interaction is taken into account. However, depending on the desired accuracy, it 
may happen that the numerical resolution, and hence computational time, required 
when using volume penalization grow disproportionately and even completely offset 
the gains. To afford a reasonable way of deciding whether one should implement a 
penalization method and, if so, how to tune its parameters, realistic order estimates 
are therefore of primary importance, taking into account not only the penalization 
error but also the discretization error. 

Even in the most simple case, which is the penalized Poisson equation, such 
estimates are either still lacking or not widely understood. Indeed, to see where the 
main issues come from, consider the following one-dimensional example: 

— w"(x) = m 2 sin ma:, (1-1) 

where the unknown w is a twice differentiable function on f2 =]0, tt[ satisfying homo- 
geneous Dirichlct boundary conditions, w(0) — w(tt) — 0, and m is an integer. The 
exact solution is 

w(x) = sinmx, (1-2) 

but forgetting it for a moment, let us approximate it using the volume penalization 
method. To this end, we extend the original domain f2 by using its natural embedding 
into the one-dimensional torus T = M./2irZ, and then look for a function on T satisfying 

— v" H — x« = m 2 sinmx, (1-3) 

where r\ > is the penalization parameter, and the mask function \ is defined by 

ifx€]0,7r[ 
x(x) = { 1/2 ifx = Oorz = tt (1.4) 

1 otherwise. 

The discontinuity of \ is the essential mathematical feature of such penalized prob- 



lems, implying in particular that the unique weak H 1 solution to ( 1.3 ) is not a classical 



solution, which in turn makes the problem hard to discretizc efficiently. Note that, 



since \ — on ]0, tt[, the restriction to ]0, w[ of any solution to ( 1.3 ) solves the original 



Poisson equation (1.1) exactly. Its discrepancy with the exact solution of the full 
Dirichlct problem is therefore entirely due to its non-zero boundary values. 

Taking advantage of the fact that the solution can be written analytically in this 
simple case (see Appendix A), one can show that the leading order L 2 error E\ with 



respect to the solution (1.2) of the Dirichlet problem is: 



mJ2 - (-1)™ ^ , . 

ei V ^ ' VV as ^0 (1.5) 

where the ^ff\ behavior is consistent with previous studies [TJ [7] . 



In order to discretize (1.3), it is possible to use a Fourier pseudo-spectral method, 



where collocation on a regular grid of TV points is used to evaluate the product \v. 
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The L 2 error between the discrete solution and the exact solution of the penalized 
problem can then be estimated using a mixture of analytical and numerical techniques 
(see Appendix A), and yields 

77i7r^/^ 1 

62 ~ K 3^2 ^ V ^ N ' (L6) 

where K = 2 for m even and K w 3.8423 for m odd. The TV -2 behavior is related to 
the regularity of the exact penalized solution, and was indeed observed for the same 
reason in a study of collocation schemes for elliptic equations with discontinuous 
coefficients [T7] , 



Combining the two estimates ( 1.5 1 and ( 1.6 1 thanks to the triangle inequality, we 
obtain a bound for the total error e between the discrete-penalized solution and the 
exact solution of the original Dirichlct boundary-value problem: 

mJl - (-l) m tott 3 / 2 1 , , N 

e<e 1 +e 2 V ^= ^ +K ^/¥^W M ? ' ~* °' ^ °°' ( ^ 

When r\ is chosen with the right order of magnitude, i.e. r\ cx 1/iV 2 , in order to 
optimize the preceding estimate, then the resulting error is 

ecxi (1.8) 

which suggests that the penalization method with Fourier discretization is, for prac- 
tical matters, a first order method, a simple fact that has been mostly overlooked by 
the literature on the subject. 

Motivated by the results obtained for the ID Poisson problem, we have sought 
to extend them and understand whether they should be seen as general properties of 
penalized equations or merely as a special case. For example, one may ask what the 
influence of the terms on the left-hand side of the equation, namely the Laplace oper- 
ator, is on the accuracy of the penalization scheme. Since one of the main domains of 
application of penalization methods is fluid dynamics, we have also oriented ourselves 
towards vector-valued penalized problems, and more specifically the penalized Stokes 
equation, which plays a role in the penalized incompressible Navier-Stokes equations. 

Another question concerns the dependency on the right hand side of the equation. 
To obtain a general answer, which characterizes the discrepancy between the original 
and penalized operators over the full range of scales, our approach in the following is 
to compare their eigenvalues and eigenf unctions. Note that the eigenvalue problem 
in itself may also be of interest for some applications. In fact, the spectral analysis 
of the penalized Laplace operator is analog to the analysis of the semi-classical limit 
for the Schrodinger operator [S] in a finite potential well, familiar from quantum 
mechanics textbooks. The penetration of the eigenfunctions within the penalized 
region corresponds to the quantum tunnel effect. By computing the convergent modes 
and monitoring their convergence rates, we obtain precise quantitative information 
on the effective accuracy of penalization methods, as well as on the properties of the 
inevitable spurious modes fostered by the extension of the domain. Then, by carrying 
out the same procedure after discretization, the properties of several schemes can be 
compared, providing realistic error estimates that can be used to evaluate the benefits 
of a penalization scheme. 

To achieve these objectives, the present work is organized as follows. First, the 
eigenvalue problems themselves are formalized, and some important notations are 
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introduced. Then the eigenvalues and eigenfunctions for the penalized Laplace op- 
erator in a ID interval are derived using a semi-analytical method. A numerical 
study of their behavior as the penalization parameter tends to zero is then conducted. 
After that, a similar analysis of the penalized Stokes operator, defined for simplic- 
ity on a straight 2D channel, is performed. A relationship between the penalized 
eigenfunctions and Navier boundary conditions is outlined. Subsequently, several dis- 
cretization schemes are compared by studying the eigenfunctions and eigenvalues of 
the discretized-penalized operators numerically. Thanks to an analysis of the penal- 
ization and discretization errors, efficient ways of choosing the parameters for realistic 
applications of the penalization method are proposed. Finally, the practical relevance 
of the preceding results is further supported by analyzing the convergence of the 
penalization method for a 2D Navier-Stokes test problem. 

2. Mathematical setting. For fi a regular bounded open subset of M. d , we shall 
denote respectively by A° and A° the Laplace and Stokes operators on fl, with their 
classical definitions: 



and, with Hp and Vq being the closures in L 2 (il) d and H 1 (Jl) d (respectively) of 



where P , called the Leray projector [TU], is the orthogonal projector on H seen as a 
sub-Hilbert space of L 2 (VL) d . 

A° and A° are positive, self-adjoint operators, with compact inverses, and there- 
fore one can construct Hilbert bases in L 2 (Q) with their eigenfunctions, which we 
denote respectively by (V>°)i£N* and (u°)i € N*- The corresponding eigenvalues, sorted 
in increasing order, are denoted by (A°)i S N* and (/i°)i S N*- 

Now assume that, thanks to a suitable similarity transformation, Q has been 
identified with a subset of the d-dimensional torus T d = (R/Z) d . The penalization 
method then consists in approximating A° by the operator A* defined on T d as follows: 



where 77 > is the penalization parameter, and \ is the indicator function of S = 
T d \Q, also called mask function. The term ^x4>i with 77 <C 1, is called penalization 
term. Note that, here and in the following, the * in the exponent serves to distinguish 
penalized operators and related quantities, while the o stands for Dirichlet boundary 
conditions. Letting 




(2.1) 



{ue Cg°(fi) I V-u = 0}, 




(2.2) 




(2.3) 



H = {ue L 2 (T d ) d I V • u = 0}, V = H n H 1 (T d ) d , 



(2.4) 



the penalized Stokes operator is defined in a similar way as follows: 




(2.5) 
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where P is the Leray projector from L 2 (T d ) d onto H. 

Like their exact counterparts A° and A", A* and j4* are linear mappings with 
compact inverses and thus enjoy Hilbert bases of eigenfunctions, denoted respectively 
(V^,j)ieN* and (u* J ieN *, while (A*^) ieN « and (M*^) iS N* are the corresponding eigen- 
values sorted in increasing order. Whereas, according to classical theorems, the eigen- 
functions of the Laplace and Stokes operators are analytic in if Q is regular enough, 
this property is not enjoyed by the penalized eigenfunctions, due to the discontinu- 
ity of the mask function. Indeed, although by construction the eigenfunctions are in 
H 2 (T d ) d , they are not in C 2 (T d ) d as we shall see below. In the upcoming sections, 
the relationship between the spectra of A° and A* on the one hand, and of A° and 
A* on the other hand, are studied using analytical and numerical techniques. 

3. Spectral analysis of penalized operators. 

3.1. Laplace operator. In the case of the Laplace operator, a ID analysis is 
sufficient to obtain a qualitative understanding of the behavior for sufficiently simple 
domains. We therefore proceed to compute the eigenfunctions of A* for d = 1 and 
f2 =]0,7r[. As a side remark, we may note the analogy with the classical eigenvalue 
problem for the Schrodinger equation (see [HI p. 67]): 

h 2 

- — Aip + ViJj = Eip, (3.1) 
2m 

where #- corresponds to ij, and the potential V to the mask function \- The limit 
77 — ^ that we study here is then analogous to the semi-classical limit h — > 0. 

3.1.1. Semi-analytic computation of eigenvalues and eigenfunctions. In 

the ID case, the eigenvalue problem reduces to a reconnection problem for ordinary 
differential equations. Indeed, the equations satisfied by ip* i and A* i can be rewritten 
as follows: 

f-r = Xtp on ]0,tt[ 

\-r= (A-i)^on ]n,2n[. [ ^ 2) 

From the form of the second equation, it appears immediately that the behavior 
will be entirely different according to whether A < rj^ 1 or A > T) . Let us first assume 
that A < 77 -1 . Applying the classical theory, we look for C 1 solutions under the form 



ip(x) 



A sm (g (x - §)) + So cos (q (x - § )) for x e ]0, w[ 
Ax (e 91 ^- 2 ^ - e^-*)) + Bx ( e q ^ x - 2 ^ + e^-*)) for x e ]ir, 2vr[ , 

(3-3) 

where q = -\/A, qx = yj r\- x — A, and the particular form of the terms has been tailored 
to avoid ill-conditionning of upcoming matrices. Since the symmetry with respect 
to a; = I leaves the problem unchanged, it should either leave the eigenfunctions 
invariant or switch their sign, which means that we can isolate two families of solutions, 
having either Bq = B\ = or Aq = A\ = 0. 

The C 1 reconnection conditions at x = n yield the following systems of equations 
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for the antisymmetric and symmetric cases, respectively: 



U sin(a)+A(l-e-^) = , . 

\A q co S (^)-A iqi (l + e-^)=0, {6 ' ' 



fl co B (^)-B 1 (l + e-^) = 
< -BbftBin(^)+fl 1 « 1 (l- e -^)=0. 



Explicit computation of the determinants of these systems yield the following eigen- 
value equations: 



F-M = sin I ^ | + Jr^-ooa [ ^ ] tanh f " /H A ] = , (3.5a) 




I-77A \ 2 / V 2 V^ 



F* + (\,r,) = cos(^ ] +^/^sinf^jcotanh^yi^A) =0, (3.5b) 
and the vanishing of F±(\, 77) is then equivalent to the existence of a non-zero C 1 



solution to ( 3.2 ), or in other words to A being an eigenvalue of A*. For fixed 77, F± (A, 77) 
are oscillatory functions of A which both admit an infinite sequence of positive roots. 
After elimination of the root A = of F* which leads to a vanishing tp and is therefore 
not an eigenvalue, the remaining strictly positive roots constitute the eigenvalues A* i 
that we were looking for. With A solution of F± (A, 77) = 0, admissible values of Aq, A\ 



or Bo,Bi are then readily obtained by solving (3.4). In the present case, this can 
be done analytically very simply, but it is not useful to give the resulting expressions 
here, since we are mostly interested in the behavior as 77 — ¥ which is studied further 
down. The final values of the constants are obtained after renormalizing the functions 
in the L 2 norm, finally yielding ip* , 



In the case A > 77 , the Ansatz (3.3) remains valid but with a complex value 



of <7i, (?i = i\J A — 77 -1 , where 1 = \J—\. Therefore, the corresponding eigenfunctions 
oscillate in ]ir, 2n[ as well as in ]0, ir[, cannot approximate the Dirichlet eigenfunctions, 
and should be considered as spurious modes generated by the penalization method. 

3.1.2. Behavior of eigenvalues and eigenfunctions as 77 — > 0. As expected, 
when 77 — > for a fixed A, the eigenvalue problem for the Laplace operator in f2 =]0, tt[ 
is recovered, with the determinants 

Fl (A) = lim F*_ (A, 77) = sin ( ^ | , (3.6a) 

F?(A) = KmnM = cos {^J > (3.6b) 

which, excluding A = as above, admit the sequence of roots A° = (i + l) 2 (i € N). 
Moreover, using symbolic computation, we have obtained that for all A > 0: 

which, by the implicit function theorem, implies that 

A;. ( = A° ( 1 - 4 ^ + 0( V )) when 77 -> 0. (3.8) 



FlG. 3.1. Eigenfunctions of the penalized Laplace operator with Q =]0, w[, for i = 1, 2, 3 (left 
to right) and varying n. Zooms close to the location of the wall are shown in the insets. 



This estimate can then be used to derive the error on the eigenfunctions. For 
this, we first establish the following useful lemma: 

Lemma 3.1. Let / :Rx!)4 R d be a continuous function which is C 1 in the 
first variable. When £ — ¥ £o; the following asymptotic expansion holds: 



/(Co, 



ll/(£,OII ||/(£o,-)ll 



l/(6,0ll ; 



(3.9) 

Proof. The relationship follows immediately from a direct computation: letting e = £, — £,0, 
we have 



n/tt,oir -«/(&, of = ey %(f)te>, •)+«(*) (3.10a) 
ii/ce.oir'ii/c&.oir^-fii/Ko.oir 3 / ^(^(fo.o+oCe) (3.10b) 

/(€,*) /(&,*) _ 9e/(Co,-)ll/(Co,-)ll 2 -|/te^)./^(/ 2 )(eo,-) 



11/&011 11/(6, 011 



i/te»,oii 3 



+ o(e) D 

(3.10c) 



Since in the antisymmetric case we have d qo ip = Aq [x — cos [q^ [x — §)), we 
obtain from Lemma 3.1 and (3.8) that 



TJ,l 



Cll = 



, 7T 1 



1 - 



(3.H) 



By evaluating |A* 4 — A°| numerically using the values of A* , provided by Newton's 
algorithm, we recover, as expected, the scaling (3.8), see Fig. 3.2 (left). We observe 
that, for low mode numbers and sufficiently small rj, the eigenfunctions tp* i of the 
penalized Laplace operator closely approximate the corresponding eigenfunctions of 
the exact Laplace operator. The cases i = 1,2,3 are presented in Fig. |3.1| In fact, it 



follows from Ansatz (3.3 ) that the eigenfunctions extend into the penalized region over 

' ' 1/2 

a boundary layer of characteristic width ^/fj (l — r/X* 4 ) . Quantitative estimates of 
the difference can be m easured in the L 2 norm by adaptive Gauss- L obat to quadrature. 
As predicted by (3.11), the error inside SI behaves as rj 1 / 2 (Fig. 3.2 middle), and 




FlG. 3.2. Left: distance between eigenvalues of the penalized and exact Laplace operators 
with Q =]0,7r[, as a function of ^/fj. Middle: L 2 distance between eigenfunctions of the penalized 
and exact Laplace operators with £1 =]0, n[ as a function of ^/rj, measured in £2 =]0,7r[ Right: same 
in S = [n, 2tt] . 



increases with the order of the eigenvalue. Inside S, the L 2 error scales like rf^ 
(Fig. |3.2| right) . Both scalings are consistent with the predictions of a boundary layer 
analysis [5]. 

The qualitative picture one can get from these results is that penalization controls 
the eigenfunctions associated to sufficiently small eigenvalues, namely those such that 
A < T) , by making them converge to the corresponding eigenfunctions of the exact 
operator with Dirichlet boundary condition. The convergence is however not uniform 



in A, since it breaks down proportionally to A°, as shown by (3.8 1 



Due to the second equation in (3.2), the remaining eigenfunctions of the penal- 
ized Laplace operator, namely those whose associated eigenvalues satisfy A > rj^ 1 , 
oscillate both in f2 and in S. They arc not related to eigenfunctions of the exact 
Laplace operator, but are spurious modes introduced by the penalization method. 
Note that (A°) has the dimension of a length and, for sufficiently simple domains, 
corresponds to a characteristic variation length of the associated eigenfunction 
Therefore, penalization allows us to control what happens over length scales coarser 
than yFq, and not finer. We come back on the issue of the choice of the right cut-off to 
eliminate these spurious modes further down in the section dedicated to discretization. 

3.2. Stokes operator. The Stokes operator makes sense for dimensions greater 
or equal to 2. For simplicity, we consider a 2D channel embedded in T 2 , £1 = ]0, ir[ x T, 
and proceed to compute the eigenfunctions of the Stokes operator in this geometry. 
For convenience, the first and second coordinates in T 2 are denoted x and y, and the 
corresponding partial derivative operators by d x and d y . 

3.2.1. Semi-analytic computation of eigenvalues and eigenfunctions. 

The method for the computation of the eigenfunctions is very similar to the one 
which we applied in the ID case above, with some additional complexities which we 
briefly describe here. The equation A*u — /iu can be rewritten: 




(3.12) 



where the effect of the Leray projector has been taken into account by introducing a 
pressure field p, which is by definition in ff 1 (T 2 ). 

Since the problem is translation invariant in the y direction, the first step is to 
expand u in terms of Fourier modes with respect to y. Denoting by u and p the 
partial Fourier tranforms of the dependent variables with respect to y, and by k the 
corresponding wavenumber, we obtain 



(d x - k z )u x - d x p = ( -[i + -x ) u x 



id 2 — k 2 )u y — ikp - 

d x u x + iku v = 0, 



V 

-A* + -X J " 



(3.13a) 

(3.13b) 
(3.13c) 



where we have used the fact that x does not depend on y. Note that the partial 
derivatives with respect to x can now be understood as ordinary derivatives, since the 
dependency on y has been eliminated. 

In the special case k — 0, (3.13b I reduces to the eigenvalue problem for the 
penalized Laplace operator in ]0,7r[, which we have already treated in the previous 
section, and (3.13c) implies that u x vanishes identically. Assuming from now on that 
k > 0, we use (3.13c) to solve for u y : 



k 



(3.14) 



and then substitute in (3.13b) to obtain the Fourier coefficients of the pressure p 

1 \ 1 



-x. 



k< 



d x u x 



(3.15) 



We now restrict our attention to either one of the two intervals ]0, tt[ and ]ir, 2tt[ over 
which x takes a constant value which we denote by k, where n = or 1 respectively. 
We then insert the expression for the pressure into (3.13a): 

1 /„, ,o 1 



d. 



H K I .' / , - 

V 



k 2 



0. 



k 2 + [i k ) d x u x = 0, 



so as to finally obtain a fourth order ordinary differential equation for u x : 



diu x 



2k 2 



1 

— / 
V 



1 



I'- 



ll 



K } U T = 0. 



By writing the solutions under the form e Lmx we obtain: 

1 

— / 
V 



in ' + j 2k 2 - fi+ -k ) iir /,-- j ir - ii 



= 0. 



(3.16) 



(3.17) 



(3.18) 



Both k and —k are roots of (3.18) whatever the value of k, and for k = 0, there are 
two other conjugate purely imaginary other roots, which we denote by tqo and — tgo • 
Now assuming that - > fi — k 2 , for n = 1 the two opposite remaining roots are real, 
and will be denoted by q\ and —qi. The values of qo and q\ are easily obtained by 
factorizing the polynomial equation (3.18): 



go = V 7* — fc 2 and qi = \J f] 1 - fj, + k 2 . 
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(3.19) 



As in the Laplace case above, the eigenfunctions have to be either symmetric or 
antisymmetric under the transformation x — > ^ — x in T, since the latter leaves \ 
invariant. Therefore we look for solutions under the form 

J A sin (g (x ~ § )) + C sinh (k (x - §)) for < x < ir 

U ' x X ' ~ [A x (e qi{x - 2n) - eM*-^) + d sinh (k (x - § )) for vr < x < 2tt, 

(3.20a) 

Bo cos (go (x — §)) + D cosh (k (x — f )) for < x < n 

Bi ( e qi(x-2TT) + e qi{-K-w)} + £) 1 cosh ( fc ^ _ |)) for 7T < Z < 2tt. 

(3.20b) 



v+{x) 



The constants are determined by the C 1 reconnection conditions for u at x = it, 
and by the continuity of p. From the expressions of u y and p as functions of u x (see 
Eqs. 3.14 and 3.15| above), these conditions may all be reformulated in terms of u x 
only: 

u x (ir + ) = u x (tt~) (3.21a) 
d x u x (n + ) = d x u x {n~) (3.21b) 
d 2 x u x {n+) = d 2 x u x (n-) (3.21c) 

d x u x {n+) = d x u x {n~) + -d x u x {TT~). (3.21d) 

V 

As before, the existence of a non-zero solution is equivalent to the vanishing of a 
determinant G± (k, /i, 77), whose expressions in the antisymmetric and symmetric cases 
(respectively) are 

G*_(k, /i, 77) = k cotanh^ — J — (1 — 277/i) (1 — 77/1) go cotan^-^J 

+ M\A(1 - ml) (1 - 2 w) cotanh (-^y 1 ) = (3.22a) 



G*+{k-,^,ri) = fctanhf^^ + (1 - 2r//j) (1 - 77/j) g tan (-|^) 



+ m^KI - Wo 2 ) (1 - 277 M ) tanh (-^) - 0. (3.22b) 

The roots of these determinants can be identified with an increasing subsequence of 
eigenvalues for A*, which we index from now on by a non-zero integer denoted I. The 
complete sequence of eigenvalues of A* n is thus indexed by the two integers k and I, or 
equivalently by i if they are sorted in increasing order as a whole. For each of these 



eigenvalues, the corresponding under-determined system (3.21) can then be solved to 
find the constants, as in the ID case. Each eigenfunction is thus computed with the 
desired accuracy. The results are shown for some values of k and I in Fig. |3.3| 



3.2.2. Behavior of eigenvalues and eigenfunctions when 7/ — > 0. When 
T] — > 0, the eigenvalue equations in the antisymmetric and symmetric cases respectively 
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FlG. 3.4. Convergence of the penalized Stokes operator in a channel A*, as function of y/rj, 
for k = 1, 2 and I = 1,2, 3, to the exact Stokes operator with Dirichlet boundary conditions A". Left: 
distance between eigenvalues. Middle: distance between eigenf unctions in L 2 (C1). Right: distance 
between eigenfunctions in L 2 (S). 



simplify to 

nk 



G°_{k,pt) = fccotanh^ — j — q cotan^-|^ = 0, (3.23a) 
G + (k,n) = fetanh^) +9otan(^) =0, (3.23b) 
whose solutions, the eigenvalues of A°, are determined numerically using Newton's 



algorithm. The corresponding eigenfunctions are computed using Ansatz (3.20) as 
above, but enforcing A\ = C\ = B\ = D\ = 0, and using the Dirichlet boundary 
conditions at x — 7r to find the values of A , Co or B , D . The eigenvalues and 
eigenfunctions of A° can then be used as reference to study the convergence of //* i 
and u*i, respectively. 
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Vv 



FlG. 3.5. \ u x ri i( n )\ (l e ft) an d \ u y r; z( n )\ (right) as functions of y/rj, for k = 1,2 and I = 1,2,3. 



As before we get by the implicit function theorem that, as r\ — > 0, 



where 



7T 



(3.24) 



1 ^cotanhC^) 



Tik\ | k A 

Mfc,,sinh 2 (^) 



1 - -^-tanh ( ^) k \, 

*7*fc,i V 2 7 M ° ,cosh 2 (^) 



for ^ even 



for I odd. 



(3.25) 



The convergence of eigenvalues thus occurs at a rate y/rj, as examplified in Fig. 3.4 
(left). Similarly to the Laplace case (Eq. 3.8), the prefactor degrades for increasing- 
mode number. To evaluate the convergence of the eigenf unctions themselves inside 
f2, the L 2 error is computed by numerical quadrature from the Ansatz (3.20) using 



the available numerical values of the parameters. The observed scalings with r\ are 
very simi lar to what we observed for the Laplace eigenf unctions, namely ry 1 / 2 in Q 
middle), and rj 3 ^ 4 in 5* (Fig. 



(Fig. 



3.4 



3.4 



right) . 



An interesting consequence of the terms C\ and D\ appearing in the Ansatz ( 3.20 1 



is that the characteristic penetration length of the penalized eigenfunctions inside the 
penalized region does not, strictly speaking, go to zero with 77, in contrast to what we 
have seen above in the case of the Laplace operator. In the case of the Stokes operator, 
the penalization error in the boundary is not limited to a vanishingly small boundary 
layer. This is due to the non-local effect of the pressure field. The effect is barely 



apparent when looking at the profiles of the ground state for varying rj (Fig. 3.3), 
because the eigenfunctions actually behave like 0{rf) deep inside the boundary. 

3.2.3. Interpretation of the leading order error as a residual slip. Since 
the eigenfunctions of the penalized Stokes problem also satisfy an eigenvalue equation 
for the Stokes operator inside f2, their discrepancy with the exact eigenfunctions of 
the latter is due to their nonzero boundary values, which we now analyze. In our 
case, the simple geometry allows us to eliminate the dependency on the wall parallel 
coordinate y simply by square-averaging in this direction. The resulting values of 

It appears that u x (ir) scales like 77 



x,r\,i 



(ir)\ and 



^71")! are shown in Fig. 
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s/V s/V s/rj 

FlG. 3.6. Left: \a* i/y/rj — l|, where a* i is the effective slip-length as a function of ^/fj for 
various eigenf unctions of the penalized Stokes operator in a channel. Middle: difference between 
eigenvalues of the penalized Stokes operator and corresponding eigenvalues of the Stokes operator in 
a channel with Navier boundary conditions and a slip length a = y/fj, as functions of n, for k = 1, 2 
and I = 1,2, 3. Right: same for eigenvalues in L 2 (Q). 



while Uy(Tr) scales like ^Jfj. In other words, the normal boundary condition of the 
limiting problem is approached much more rapidly in the penalized problem than the 
tangential boundary condition. This observation, as well as some observations that 
we made when solving the Navier-Stokes equations [2U] , has prompted us to introduce 
the Stokes eigenvalue problem with Navier boundary conditions, as we now present. 

The so-called Navier boundary conditions, already introduced by Navier in the 
historical paper where the equations later known as the Navier-Stokes equations were 
proposed [19], can be written 

uan • t + ■ r = (3.26a) 

on 

u dn • n = 0, (3.26b) 

where n is the outward normal unit vector, r is any tangent unit vector, both along 
<9f2, and is a positive real number known as the slip length. Dirichlet boundary 
conditions, which we have been considering up to now, correspond to the special case 
cr = (no slip). Coming back to the penalized problem, we can thus define, by 



analogy with (|3.26a|), an effective slip length a* v for its eigenfunctions: 



assuming that the denominator does not vanish. The surprising fact that we have 
found is that a* „ is close to y/rj independently of i. This can be seen for example by 
monitoring the difference {a* irj — v^Vy^, which is found to scale like 0(y/rj) (Fig. 



3i)| left). 

This finding suggests a direct comparison between the eigenfunctions of the pe- 
nalized Stokes operator on the one hand, and those of the Stokes operator inside Q 
supplemented with Navier boundary conditions on the other hand. For this, define 
Vt, the closure in iJ 1 (fi) rf of 

dxx 

u € C°°(0) I V • u = 0, u an • n = 0, u an • t + • r = 

an 
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and consider the operator : (H 2 (fi)) nV ' v — > H, which is defined exactly as A° but 
by replacing Vo by in the definition of the starting space. The eigenvalue problem 
for is then well defined, and the corresponding eigenvalues and eigenfunctions will 

be denoted respectively by fi\ and uj . They can be computed semi-analytically in 
the same manner as was done for A° above, simply by enforcing the Navier conditions 
(3.26) at the boundary instead of Dirichlet conditions, thus yielding both different 
eigenvalues and different eigenfunctions. The corresponding eigenvalue equations are: 



Gl(fc, fi,r]) = fccotanh^^-J — q cotan^^-^ — ^/rjfi = 0, (3.28a) 
Gl(k,fi,rj) = fctanh^j + q tan(^) + ^ = 0, (3.28b) 



to be compared with (3.23) and (3.22) 



As we conjectured, the numerical evaluation of \[i* v — fJ-l v \ as a function of ij 



suggests that it is proportionnal to 77 (Fig. 3.6 middl e), a nd the distance between 
the eigenfunctions in L 2 (fl) is of the same order (Fig. 3.6 right). This supports a 
leading order description of the error associated to the penalized eigenfunctions by 
the error coming from the residual slip. This description has important consequences 
for applications of the penalization method in computational fluid dynamics, and it 
would be of interest to derive a proof for general geometries. In the present case, the 
leading order estimate obtained by the implicit function theorem is the same as for 



the penalized eigenvalues (Eq. 3.24), which proves that 



Av = ^tn + ^)- ( 3 - 29 ) 

4. Discretization. In this section we compare several numerical methods that 
can be used to discretize the eigenvalue problems for the penalized Laplace and Stokes 
operators. 

4.1. Laplace operator. 

4.1.1. Numerical methods. The discretization of penalized operators is not 
an easy problem because they involve discontinuous mask functions. In particular, 
when considering the eigenvalue problem in Fourier space, a way must be found to 
truncate the Fourier series of the mask function, or regularize it somehow. Here, we 
compare the following approaches: 

(i) a Fourier collocation method, where the product between the (discontinuous) 
mask and the function is evaluated over a regular grid in physical space, 

(ii) a Fourier-Galcrkin method, with a sharp truncation of the discrete Fourier 
series of the mask function similar to [35] , 

(iii) a Fourier- Galerkin method, with a smooth, positivity-preserving truncation 
of the discrete Fourier series of the mask function (as detailed in Appendix[B|), 

(iv) second or fourth order centered finite differences. 

Note that the convergence of the collocation scheme (i) for an equivalent eigen- 
value problem (but only for fixed rf) was already studied in |17j . In order to understand 
the convergence rate of the collocation scheme, the authors proved estimates concern- 
ing an "ideal" Fourier-Galerkin scheme, using the exact Fourier coefficients of the 
mask function, and obtained a convergence rate of —2.5. This should not be confused 
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FlG. 4.1. L 2 distance between the first eigenf unctions of the discrete-penalized Laplace operator 
and of the penalized Laplace operator, as a function of N and for varying n. Top: Fourier collocation. 
Middle: Fourier- Galerkin with sharp truncation. Right: Fourier- Galerkin with positive mollified 
mask. 



with the problems (ii) and (iii) that are under consideration here, which rely only on 
the knowledge of the discrete Fourier coefficients of the mask function. 

The Fourier-Galerkin methods are implemented by the pseudo-spectral approach, 
which consists in evaluating the product x^V' i n physical space. To allow a computa- 
tion of the Fourier coefficients of without any aliasing error, we take a grid having 
N = 4K points in each direction, where K is the cut-off wavenumber. The FFTs are 
performed using the FFTW library [T2]. Whatever the discretization method, we 
obtain a finite dimensional eigenvalue problem, which we then solve numerically by a 
Krylov iteration method, as efficiently implemented in the Petsc and Slepc libraries 

In the following, we present the results of a parametric study where 81 values of 
77, logarithmically equidistant between 10 -2 and 10 -6 , and 16 values of N, of the form 
N = 2% N = 3 x 2 J_1 where 6 < i < 13, have been screened. 

4.1.2. Results. First, we examine the relationship between the discretized- 
penalized Laplace operator on the one hand, and the continuous penalized Laplace 
operator on the other hand. In the collocation case, the first eigenfunction converges at 



a rate N 2 (Fig. 4.1 left), consistent with the results proved in [T7]. This convergence 
rate is connected to the regularity of the limiting eigenfunction (see also Appendix [A]) 



4.1 



The same convergence rate is observed in the sharp Galerkin truncation case (Fig 
left). Hence the optimal TV -2 5 convergence rate which would correspond to the ideal 
Fourier-Galerkin scheme without truncation of the mask function, and using its exact 
Fourier coefficients (see [T7]), is not attained here. We have checked that the same 
behavior is observed for higher order eigcnfunctions (up to i = 5). In the smooth 
Galerkin-truncation case, the mollification of the mask introduces perturbations of 
order A^ 1 in the equation, which results in perturbations of the same order on the 



eigenfunctions (Fig. 4.1 right). The error displays a faster decay for sufficiently small 
values of N. 

We now turn to the distances between the first eigenfunction of the discretized- 
penalized Laplace operator and of the exact Laplace operator with Dirichlet b.c, as 



depicted in Fig. 4.2 as a function of ^Jfj. As expected, for small r\ the 0(y/rj) behavior 
of the penalization error dominates. Interestingly, for all three schemes the error 
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FlG. 4.2. L 2 distance between the first eigenf unctions of the discrete-penalized Laplace operator 
and of the exact Laplace operator with Dirichlet b.c, as a function of Jf) and for varying N . For 
comparison, the dash-dotted line recalls the error for the exact penalized eigenfunction. Left: Fourier 
collocation. Middle: Fourier- Galerkin with sharp truncation. Right: Fourier- Galerkin with positive 
mollified mask. 



locally drops down to a minimum which is of order ?y 3//4 . Then, for the collocation 
scheme a plateau is observed, where the error saturates, whereas for the two other 
methods, the error seems to increase without bound for smaller and smaller rj. The 
presence of the rj 3 ^ 4 minimum came as a surprise to us, since it cannot be explained 



by naive bounds such as (1.7) given by the triangle inequality, and must be due 
to unexpected cancellations between the penalization and discretization errors. We 
comment on its significance further down. 

In Fig. |4.3| (left) the error associated to all considered schemes is plotted as a 
function of ^fr) for a fixed value of N, which allows for a direct and quantitative 
comparison. In addition, the errors obtained when using 2nd and 4th order finite 
differences schemes are shown. The collocation, truncated-Galerkin, and 4th order 
finite differences schemes have similar properties, while the smoothed-Galerkin scheme 
saturates at a value of the error which is roughly 10 times larger. The behavior of 
the 2nd order finite difference scheme stands out, since the range of rj where the error 
saturates is not even reached in our study. 

In practice, rj is an adjustable parameter which should be chosen in order to 
minimize the error with the eigenfunction one wishes to approximate. We denote by 
Vopti such a well-chosen value of rj. Its behavior as a function of N, and the actual 



value of the error attained for rj = rj op ti, are shown in Fig. 4.3 (middle) and Fig. 4.3 
(right ) , resp ect i vely. 

It appears that rj op ti scales like N~ 2 , similarly to what we found for the Poisson 
problem presented in the introduction. However, the resulting error when choosing 
such an optimal rj scales like N~ 3 / 2 and not iV _1 as we had expected from the Poisson 
example. Note that the range of t) for which these cancellations occurs is relatively 
narrow and may be difficult to estimate in many cases of practical interest. 

4.2. Stokes operator. 

4.2.1. Numerical method. In the case of the Stokes operator, we leave aside 
finite differences due to technical constraints, and we concentrate on the collocation, 
truncated-Galerkin and smoothed-Galerkin approaches. The problem is first reduced 
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FlG. 4.3. Comparison of various discretization schemes. Left: L 2 distance between the 
first eigenf unctions of the discrete-penalized Laplace operator and of the exact Laplace operator 
with Dirichlet b.c, as a function of ^/rj, for N = 256. Middle: Value rj opt i of r\ at which the 
minimum error between the first eigenfunctions of the discrete-penalized and exact Laplace operators 
is attained, as a function of N . Right: corresponding minimal L 2 error. 



to a scalar problem by taking the curl of the eigenvalue equation, which yields an 
equation for the vorticity uj of the eigenfunction: 

-Acj + - (d x {x*u y ) - d y {x*u x )) = [MJ 

The dependency on y is then eliminated by taking the Fourier transform in the y 
direction, after which all modes with different wavenumbers become decoupled, and 
k becomes simply a parameter. We thus obtain a set of uncoupled one-dimensional 
problems which, although not Hermitian symmetric, can be solved for up to N = 2 12 
on a standard personnal computer within a few hours using, as above, the Petsc/Slcpc 
libraries. Once this is done, the velocity field is reconstructed from u> by applying the 
Biot-Savart operator in Fourier space. Otherwise, the procedure and the values of ./V 
and 77 are the same as in the Laplace case as reported above. 



4.2.2. Results. Concerning the discretization error (Fig. 4.4), as well as the 



total error (Fig. 4.5), we do not observe any notable difference with the penalized 
Laplace operator. We have checked that this also holds for higher order eigenfunctions 
(up to k = 3 and 1 = 5). 

The comparison with the first eigenfunction of the exact Stokes operator with 



Navier b.c. with slip length ^/rj is depicted in Fig. 4.6 The error scales like 0{rf) for 
large 77, as expected from the analysis of the penalized eigenfunctions in the previous 
section, and does not have such a sharp minimum as in the Dirichlet case. 

In the case of Dirichlet b.c, the behaviors of r] opti and of the minimal error 



(Fig. 4.7 1 are analogous to what we have seen above in the case of the Laplace operator, 
namely r\ op n scaling like iV -2 and the optimal error like iV~ 3 / 2 . For Navier boundary 
conditions, we find the much weaker scaling iV~ 4 / 3 for r] opt i, which means that the 
optimal error also scales like N~ 4 / 3 . Although penalized eigenfunctions for a given 
r\ are much closer to Navier eigenfunctions with slip length ^Jf\ than to Dirichlet 
eigenfunctions, this advantage seems to be destroyed by the discretization when rj 
gets smaller than iV -4 / 3 . 
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FIG. 4.4. L 2 distance between the first eigenf unctions of the discrete-penalized Stokes operator 
and of the penalized Stokes operator, as a function of N and for varying n. Left: Fourier collocation. 
Middle: Fourier- Galerkin with sharp truncation. Right: Fourier- Galerkin with positive mollified 
mask. 
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FlG. 4.5. L 2 distance between the first eigenf unctions of the discrete-penalized Stokes operator 
and of the exact Stokes operator with Dirichlet b.c, as a function of ^/rj and for varying N. For 
comparison, the dash-dotted line recalls the error for the exact penalized eigenfunction. Left: Fourier 
collocation. Middle: Fourier- Galerkin with sharp truncation. Right: Fourier- Galerkin with positive 
mollified mask. 



The fact that the error observed with respect to Navier b.c. is higher than with 
respect to Dirichlet b.c. is somewhat disappointing. It is related to the occurence 
of unexpected cancellations in the Dirichlet case, which do not occur in the Navier 
case. A possibility which remains, to be investigated in the general case, would be 
that the discrete eigenfunctions are better described by Navier boundary conditions 
with a slip length different from ^/rj. 



5. Application: dipole-wall collision. In this section we briefly illustrate 
the connection between the behavior described above for the linear Stokes eigenvalue 
problem, and the initial-boundary value problem for the incompressible Navier-Stokes 
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FlG. 4.6. L? distance between the first eigenf unctions of the discrete-penalized Stokes operator 
and of the exact Stokes operator with Navier b.c, as a function of ^/fj and for varying N . Left: 
Fourier collocation. Middle: Fourier- Galerkin with sharp truncation. Right: Fourier- Galerkin with 
positive mollified mask. 




FlG. 4.7. Comparison of various discretization schemes. Left: value r\ ov ti of n at which the 
minimum error between the first eigenfunctions of the discrete-penalized and exact Stokes operators 
is attained, as a function of N. Right: corresponding minimal L 2 error. The solid lines corresponds 
to Dirichlet boundary conditions, and the dash-dotted lines correspond to Navier boundary conditions 
with a slip-length ^/rj. 



equations: 



d t u + (u • V)u = - Vp + vAu 
(NSE) : I V • u = 

k U| 9n = 0, u(0,x) = u o , 



(5.1) 



where v is the kinematic viscosity, and the spatial domain f2 =]0,7r[xT is as before. 
The penalized version of (5.1) is commonly written as: 



<9 t u + (u • V)u = — Vp + i/Au — 7V xu 
(PNSE) : \ V • u = 

k u(0,x) = uo, 
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(5.2) 





N 


rf 


i] := 


urj' 




Vv 




I 


1024 


0.63 


9.87- 


io- 


5 


9.93 • 10- 


3 


II 


2048 


0.25 


3.95- 


io- 


5 


6.28 • 10- 


3 


III 


4096 


0.10 


1.58- 


io- 


5 


3.97-10" 


3 


IV 


8192 


0.040 


6.32- 


io- 


6 


2.51 • 10- 


3 



Table 5.1 

Parameters used for the four considered Navier-Stokes test cases. 



where x now varies in T 2 . Note that rf in this equation has the dimension of an 
inverse time, whereas rj, the penalization parameter as we have considered it up to 
now, is a squared length. In order to connect the two quantities, we have to rewrite 



(5.2) as follows: 

d t u + P[(u • V)u] = uA*u, (5.3) 

where r\ := rj v is now the appropriately scaled parameter for the penalized Stokes 
operator as defined by (2.5), and P is the Leray projector on T 2 . 



For spatial discretization, we use the Fourier-Galerkin scheme with positively 
smoothed mask function, with a pseudo-spectral implementation on an N x N grid 
for efficiency. Time discretization is achieved with a third order low storage Runge- 
Kutta method [2TJ page 20], and the timestep is adjusted dynamically to satisfy 
the CFL condition. To analyze the influence of the penalization and discretization 
parameters on the solution of (PNSE), we shall vary rj and N only and keep all the 



other parameters fixed. N is varied by powers of 2 from 1024 to 8192 (see Table 5.1 ) 



and the values o f rj a re chosen according to the expected behavior of the optimal 77 



according to Fig. 4.7 i.e. with the scaling rj oc N~ 4 ^ 3 . 

As initial condition, we choose a quadrupole centered around x — 7r/2, y — tt, for 
which the vorticity field is given by: 

Uo(%,y) ■= d x u ^ y - d y u 0<x := -^xy{x 2 + y 2 ~ 6s 2 )e~ 1^ , (5.4) 



where A ~ 0.6258473 and s = 0.2, as shown in Fig. 5.1 (top, left). Its initial enstrophy 



is Z = 1 f n lu 2 ~ 7.38309-10" 2 , and its initial energy is E = \ J Q u§ ~ 1.87016-10" 5 . 
Oriented contours of the initial stream function ipQ, define by V 2 i/'o — w 0j ar e shown 



in Fig. [5TJ (top, left). 

As viscosity we take v ~ 1.57914 • 10~ 4 , which corresponds to an initial Reynolds 
number based on root mean square velocity and channel width of Re = \/2Eqit jv ~ 
122. Note that this Reynolds number is taken sufficiently low so that the discretization 
parameter can be varied over a sufficiently wide range while still properly resolving 
the solution. 

Vorticity contour plots are shown in Fig. |5.1[ As expected, the boundary layer 
becomes sharper as N increases and 77 diminishes. Note that, thanks to the smoothing 
of the mask function, spurious oscillations do not arise. 



Now turning to more quantitative diagnostics, we compare in Fig. 5.2 (top left) the 



residual tangential and normal velocities, in the spirit of Fig. 3.5 The same scalings 
are obtained as in the linear Stokes case, namely ||ti x |an|| oc rj for the normal velocity 
and Hu^ianll °c y/rj for the tangential velocity. By plotting the residual tangential 



velocity at the wall u v \qq as a function of the wall strain rate d x u y \QQ (Fig- 5.2 top 
right) we observe a linear correlation, as expected if Navier boundary conditions are 
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FlG. 5.1. Dipole-wall collision at Re = 122. Top left: contour plot of the initial stream 
function (which is identical in cases I-IV), with 12 contours equally spaced between —0.01 and 0.01. 
Top right: contour plot of vorticity at t = 50 for r\' = 0.25 and N = 2048 (case II), restricted to the 
subdomain [37r/4,7r]x [3tt/4,tt], with the following contours lines: —0.28, —0.21, —0.14, —0.07 (dash- 
dotted lines), and 0.02, 0.04, 0.06 (full lines). Bottom left: same but for r\' = 0.1 and N = 4096 
(case III). Bottom right: same but for r\' = 0.04 and N = 8192 (case IV). 



approximately satisfied (see (3.26a)). The effective slip-length a*, as obtained from a 

2 (bottom left). Recall 
Here, this relationship 



5.2 



least-squares fit, is then plotted as a function of ^Jrj in Fig. 
that we have shown in the linear Stokes case that a* ~ ^frj. 
seems to be approached in the limit 77 — 0. 

All these diagnostics suggest a very good agreement between the picture we have 
got from the linear case on the one hand, and the Navier-Stokes test case on the 
other hand. The limits of this approach are however seen when measuring the higher 
order residual + a*d x u y \Qci\\ + \\u x \qq\\, which is an indicator of the precision 

with which the solution matches Navier boundary conditions with the self-consistently 
determined slip length. The curve (Fig. 



5.2 



bottom right) suggests that this residual 
saturates for small values of the penalization parameter. We attribute this saturation 
to nonlinear effects. Note however that, despite this saturation, the residual remains 
about one order of magnitude smaller than the tangential velocity (Fig. 5.2 top 
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FlG. 5.2. j4 posteriori analysis of boundary conditions for the penalized Navier-Stokes test 
cases at t = 50. Top left: wall-tangential and wall-normal residual velocities as functions of ^/rj. 
Top right: wall-parallel residual velocity as a function of wall strain rate. Bottom left: effective 
residual slip-length a* as a function of yfr\. Bottom right: L 2 discrepancy along the boundary with 
respect to Navier boundary conditions with slip-length a* . 



left), suggesting that the penalized solution indeed approaches the Navier boundary 
condition faster than the no-slip one. 

6. Summary and conclusion. Through a detailed analysis of the spectral de- 
composition the penalized Laplace and Stokes operators in simplified geometries, the 
following properties were observed. 

The O^ 1 ' 2 ) bound on the error with respect to the Dirichlet boundary condi- 
tions, which was expected from previous studies, indeed applies to eigenfunctions of 
sufficiently low order. However, the convergence rate also depends on the eigenvalue, 
which entrains a degradation for higher order eigenfunctions, and even a break down 
of convergence when the eigenvalue exceeds a certain value (0{ri~ r ) in the Laplace 
case and 0(t7 —1 — k 2 ) in the Stokes case, respectively). On dimensional grounds, 
this means that the penalization does not control the behavior of the solution below 
a certain scale (O^ 1 / 2 ) in the Laplace case and 0((?7 _1 — fc~ 2 ) -1 / 2 ) in the Stokes 
case, respectively). In the case of the Stokes operator, the error can be described at 
lowest order by Navier boundary conditions with a residual slip-length proportionnal 
to O^ 1 / 2 ), which may have important consequences for the physical interpretation 
of results obtained with the volume penalization method. 

The discretization of the penalized Laplace and Stokes operators using spectral 

22 



schemes brings some new important elements into the story. In both cases, the optimal 
penalization parameter r\ ov ti is of order N~ 2 , in agreement with the notion that the 
finest resolved scale should be of order r\ x l 2 . Thanks to unexpected cancellations 
between the penalization and discretization errors, the total error with respect to the 
Dirichlet solution for r\ = r\ ov u is of order N~ 3 / 2 . These cancellations occur only for a 
narrow range of values of r\ very close to rj opt i, and we were able to observe them only 
by systematically scanning for the right value of 77. This is usually not possible in 
practical applications of the method, when the exact solution is unknown. However, 
we expect the scaling rj opt i oc N~ 2 to be very general, and we propose it as a guideline 
for future applications of the penalization method with Fourier-type discretizations. 

If one wishes to reach a higher precision with a penalization scheme, one possi- 
bility is to use a discretization which behaves better in the presence of discontinuous 
solutions. For example, a second order finite volume scheme was considered in |23j . 
involving a careful correction in the penalization term. 



Appendix A. Analytic solution and discretization of the ID penalized 
Poisson problem. 

In this appendix, we provide a detailed analysis of the introductory example 
concerning the ID penalized Poisson equation, as well as of its discretization. 

A.l. Analytic solution of the penalized Poisson equation. The exact so- 
lution to the penalized Poisson problem (1.3) is: 

{sin ma; + A\X + A 2 , x€[0,7r[ 
T fg T sinmx + B ie -^Vv + B 2 e x '^ 1 x e [n, 2ir[, 



(A.l) 



where 
Ax = 

K -- 



1 + r\m 2 
2 



K. 



A, = 



1 + r]m z 

e 2n /^[A 2 {l - e-^/v^) + Am] 



(K+l)cosh^-K-(-l) 



sinh ^= 



2 sinh 

VvQ- + (-1)" 

7rcotanIi2^ + 2y/rj 



Bo 



2 sinh J= 



(A.2) 

It is visualized in figure A.l for m — 2. The second derivative v" is discontinuous at 
x = and x = ir. The penalization boundary layer inside of the solid becomes thinner 
as i] — > and, in the limit, the first derivative v' is discontinuous. 

Knowing (1.2) and (A.l), we directly calculate the penalzation error, 



2 IT 



1 

2^ 



(v - w e ) 2 dx 



l-L 2 (]0,7T[) 



lL 2 (]7r,27r[) 



(A.3) 



where w e is the exact solution (1.2) extended to T, 

e/ . { smmx, xG [0, 7r[ 
w (x) = < „ ' , 

w \ 0, x G [7r,27r[. 
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(A.4) 




The error norm in the fluid domain is 

7T 

1 



m 2 r] 



^Q0M) = Z l(Aix + A 2 ydx~^(2-(-ir) as V -+ 0, (A.5) 



and in the solid domain 



2tt 



\L 2 (]tv,2tt[) 



1 + r\nv 



■smmx + B ie - x/ ^ + B 2 e x/ ^ ¥ > ) dx 



mW /2 



as r\ —> 0. 



(A.6) 



The penalization error in the fluid decays slower then in the solid, hence 



e 2 



y/6 



■y/rj as r/ — ^ 0. 



(A.7) 



A. 2. Discretization error of the pseudo-spectral solution. The Fourier 
coefficients of (A.l) are given by 



Vk 



2tt 



v{x)e~ lkx dx 



(A.. 



where k G Z is the wavenumbcr. This can be integrated in terms of elementary 
functions, but herein we only need (A. 8 1 expanded in negative powers of k for large 
k, and then the coefficients of that series are expanded in 7?, giving 



Vk ' 

rn 
2^ 



2tt 



-d + (-1)^)4= + (1 + ( ' 1)m)(1 + ( - 1)fc) + om 



1 - (-l) m + fe 



1 



+ (i + (-ir)(i-(-i) fe ) 



i 



i 



i 



(A.9) 
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FlG. A. 2. Absolute value of the Fourier coefficients for m = 2 and two values of the penaliza- 
tion parameter: r\ = 10~ 2 (left) and r\ = 10~ 4 (right). The solid line depicts the exact values \A. 
Gray open circles show the leading-order term in ^A.9\ . Black dots correspond to {A.llty . 



It is also instructive to consider the Fourier coefficients of (A.4|, 
1 m(l - (-l) m+fc ) 



2tt m 2 — k 2 

-t/4, 

t/4, 



k 7^ m, k =^ —m 
k = m 

k = — 771. 



(A.10) 



The leading-order asymptotic approximation to (A. 10 1 for large k is 

771 (1 - (-l) m+fc ) 



It'l 



2tt 



k 2 



+ o(k- 2 ), 



(A.ll) 



which shows that decays slower than or equivalently that w e is less regular 
than v. Fig. A. 2 shows the Fourier coefficients calculated with two different values of 
7]. The wavenumber corresponding to the slope change is estimated from (A. 9) and 
lATil as 



h = i/Vv- 

The discrete Fourier transform of v(x) is defined as 
N-x 



(A.12) 



i k = — v(x n )e~ Lkx ", where x n = 2irn/N, k = -N/2. ..N/2 - 1 (A.13) 

n=0 

v(x) is real- valued, therefore the discrete Fourier transform is conjugate symmetric, 
hence Vk = v*_ k . The inverse of (A.13) is 



JV/2-1 

v(x n ) = ^ v k e" 

k=~N/2 



(A.14) 



AT/2-1 

Note that the truncated Fourier series v k e LkXn is not interpolating, but it min- 

k=-N/2 

imizes the L 2 norm of the error with respect to v{x). Also, if the Fourier coefficients 
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Vk, as defined by (A.8), are known, the discrete Fourier transform can be obtained by 
periodization: 



Vk 



Vk+Np- 



(A.15) 



Let us now solve the penalized Poisson equation (1.3) using a pseudo-spectral 
Fourier method. The spatial domain is discretized with N grid points {x n = 2Trn/N}, 
n = Q...N — 1. The second term on the left-hand side of (1.3) is evaluated at these 
points, then the discrete Fourier transform is performed. Finally, the pseudo-spectral 
solution u is obtained from a linear system of equations for Uk 



k 2 u k + - {xu) k = fu- 
rl 



(A.16) 



We are now interested in estimating the error between the pseudo-spectral solution 
u(x) and the exact solution v(x) of the penalized Poisson equation. Let us rewrite 
\k.l&\ as 



k 2 (u k -Vk + v k ) + -({xu) k ~ (xv)k + (xv) k ) = fk- 



(A.17) 



Note that 

efe = u k - v k (A. 18) 

are the coefficients of the error between the pseudo-spectral and the exact solution, 
i.e. the discretization error. Using (1.3), the right-hand side can be expressed in 
terms of the exact solution v, 



h = ~(v") k + \ X v)k- 



From (A.17), (A. 18), (A. 19) follows the equation for the discretization error, 

k 2 e k + ^(xe) fe = -{v") k ~ ^"fe- 



lt is useful to rewrite (A.20) as 



1 



TTW [e ~ fc " {xe)k] = IT^ hK)fc " k ' 2ik 



Our objective is to estimate ||e||2 without solving (A. 21). We note that 

||e|| 2 =i^- 1/4 |H| 2 



(A.19) 



(A.20) 



(A.21) 



(A.22) 



where b is the right hand side of (A.21), and K = 2 for m even and K w 3.8423 
for m odd. This relationship was confirmed numerically with arbitrary values of the 
parameters, but it is unclear how to prove it analytically. 

Now let us derive the leading-order asymptotic approximation for 



N/2-1 

E 

k=-N/2 



V 



1 + r]k 2 
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(A.23) 



Periodization gives 



-(«")*= E {k + Np) 2 v k+Np , 

P — — OD 

OO 

-k 2 v k = —k 2 ^ Vk+N P - 



(A.24) 



We assume that N is even, hence 1 ± (— l) k + N P = 1 ± (-l) k . Using (A. 9) we obtain 

m 2 7T 2 (1 + (-ly^TJ^k 2 + (1 - (-l)"^^- 2 



18 



-(«")* 

Q( v -^ 2 )k 2 + Q(7T 3 / 2 ) 

g 0£p(ta-l))O(fc2) + O( V - 2a )O(k ) 



N 4 



a=2 
J/2 J/2 

EE 

0=1 7=1 



-(£+7-1/2) 



(A.25) 



)0(fe 1 ) 



AT2/3+27 



^ ^ O(7 7 -(^- 1 ))O(fc 2 )+O(r,-(' 3 +^)O(A ; ) 



_/V 2 / 3 + 2 7 



We notice that expansion (A. 9 1 is only valid for large k compared to rp 1 / 2 . However, 
the truncation and aliasing errors due to the spectral discretization of any function 
are determined by its high-wavenumber residual. The low-wavenumber modes Vk, 
though entering equation (A.23), must cancel out. Therefore the error committed 
in the lower part of the spectrum in (A.9) does not enter the large- N asymptotic 
estimates ( |A.27[ )-( [A~28| belo w. 

After multiplying (A.25 1 by (fc 2 + l/r))~ 2 , it is convenient to rewrite its first term 



as 



^ 2 tt 2 (1 + (-l)"**)^- 1 * 2 + (1 - (-l) m+fc )ry- 



18iV 4 



18iV 4 



1 



r\ \k 2 + 1/rj 



(fc 2 + l/ry) 2 
1 (-D fe 



k 2 + 1/rj 



2(-l) m (-1)* 



rj 2 (k 2 + l/rjy 



(A.26) 



From (JAT23J), ( |A.25| > and (JAT26J), summing up for k = -N/2. ..N/2 - 1, we obtain 

1 



/ m 2 7r 3 



+ 0(77 



1 



_/VP+3 



(A.27) 



and (A. 22) yields the desired estimate: 
.m7T 3 / 2 1 



l|e|b ~ ( A + o( " 0) ) + + ■•■ + ° ( "" /S) ^ « 6 R 

(A.28) 

Figure A. 3 presents a comparison of the first term in (A.28) with the results of the 
numerical pseudo-spectral solution of (1.3). 
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m=2, r\ =0.01 



m =2, y\ =0.0001 




FlG. A. 3. Absolute value of the Fourier coefficients for m = 2 and two values of the penaliza- 
tion parameter: r\ = 10 — 2 (left) and rj = 10 -4 (right). The line with circular markers is obtained 
from a numerical solution discussed in section \A.3\ Crosses show the first term of the asymptotic 
estimate ![A.28\) . The discrepancy between this estimate and the numerical solution is shown with a 
solid line. 



A. 3. Numerical tests. 



Convergence tests have been performed for three differ- 
a pseudo-spectral Fourier, second- 



ent discretization schemes employed to solve (1.3) 
and fourth-order finite difference methods. For the Fourier spectral method, the sec 
ond derivative operator is approximated as 

\ 2 



D 2 
U F 



/ 





— |cotan^ 




— |cotan^ 


|cotan^ 




^cotan 2 ^ 


— icotan^ 




— |cotan^ 








^cotan^ 


V 


^cotan^ 






the second-order central finite difference operator is 



D 2 - — 



(-2 1 
1 -2 1 



1 \ 



V 1 !-->/ 

and the fourth-order central finite difference operator is 

/ 



n 2 

LJ FD4 



1 



V 



_2 
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4 
3 
1 
12 



12 



12 
4 
3 
_5 
2 



_1_ 
12 

4 _X_ 
3 12 



12 



(A.29) 



(A.30) 



3 » 
J_ 

12 



(A.31) 



where h = 2n/N. The right-hand side f(x) and the mask function x( x ) are evaluated 
at the grid points x n = 2i:n/N, n = 0, ...,N — 1. The numerical solution of the 
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penalized equation is then obtained by solving the linear system 



(-D 2 + - X l)u = f, (A.32) 
V 

where I is the identity matrix and D 2 is the appropriate differentiation matrix. Com- 
parison of this numerical solution with the analytical solutions w and v is presented 
in Figs. |A.4| |A.5| and |A.6| The error e w with respect to the solution of the original 
boundary value problem w is calculated only in the fluid domain [0, 7r], and the error 
e with respect to the penalized solution v is calculated in the periodic domain [0, 2tt[. 

In all cases, two regimes are observed, depending on whether the grid is sufficiently 
dense to approximate the solution in the penalization boundary layer inside the solid. 



For the pseudo-spectral method (figure A. 4 1, the error e w decays like A 1 until N 



N t — 2/y/fj, when the spectrum is truncated at k = k t (see (A. 12 1 and Fig. A. 2 1. A 
faster convergence is observed when the number of grid points is near Nt, but with 
N increasing further on, the error starts growing and asymptotically behaves like the 
penalization error e = \v — w e \. 

The error e also exhibits two distinct regimes, but there is no super-convergence 
in between. When N is small, the observations suggest a power law e « r] a N~ 3 / 4 
with a w 0.071. When N is large enough, e decays like A -2 but diverges with 77 — >• 



like l/^frj, in agreement with (A. 3). 

Interestingly, the second-order finite difference enjoys a faster convergence in the 
intermediate regime (figure A.5). Both e w and e decay like N~ 2 when N is small. 
For large AT, e decays like N~' z and e w saturates. 

The rapid convergence in this case can be explained as follows. The second- 
order finite-difference discretization results in linear algebraic equations which are 
exactly the same for the penalized problem as for the original Dirichlct boundary- 
value problem at all points inside the fluid domain, except x\ and xj\r/2— 1> At these 
two points, an equivalent near-boundary scheme for a non-homogeneous Dirichlet 
boundary- value problem can be written as 

u 2 -2ui+a /3 - 2u N j 2 -\ + Ujv/2-2 _ , / \ oo\ 
p - h ' p - J N/2-1 , (A.66) 

where a = f3 = for the original homogeneous Dirichlet problem and max(|a|, < 
f{h)y/rj for the penalized problem. Therefore, for a fixed h and sufficiently small "q, 
the departure from the homogeneous Dirichlet boundary condition is negligible and 
the rate of convergence of 0(h 2 ) is observed in a wide range of h. The fourth-order 



finite-difference scheme performs very much like the spectral method (Fig. A. 6 ) 



A. 4. Concluding remarks. A Fourier pseudo-spectral discretization of the 
one-dimensional penalized Poisson equation has been studied. For a fixed 77, second- 
order convergence with N is found with respect to the exact solution of the penalized 
problem, as expected from its regularity. Moreover, the analysis shows that the pre- 



factor of this asymptote diverges like l/y^ as 77 — > 0, see (A.28). 

A matter of more practical interest is the convergence of the numerical solution 
towards the exact solution of the Dirichlet boundary value problem. For a fixed 77, the 
total error (of penalization and truncation) has a minimum at N = 2/- s /fj. Varying 77 
with N like 77 = 4/A^ 2 would allow the error to decay approximately like A -1 . 

Interestingly, the total error of the second-order finite-difference discretization de- 
cays like A~ 2 in a wide range of A, if there is a grid point coinciding with the bound- 
ary and the mask function at that point is greater than 0. This happens because the 
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FlG. A. 4. Convergence of the Fourier pseudo-spectral scheme, m = 2. Error with respect 
to the exact Dirichlet solution in the interior of the fluid domain (left) and with respect to the the 
penalized solution in the whole domain (right). 




FlG. A. 5. Convergence of the second- order finite difference scheme, m = 2. Error with respect 
to the exact Dirichlet solution in the interior of the fluid domain (left) and with respect to the the 
penalized solution in the whole domain (right). 



second-order finite-difference discretization of the penalized problem becomes equiv- 
alent to a non-homogeneous Dirichlet boundary value problem with negligibly small 
values at the boundaries. 

Appendix B. Positive-definite mollified penalized operators. 

Here we explain the details of the smoothing approach which was used in Section 4. 
The idea is to approximate % by a smooth function prior to discretization using a 
Fourier-Galerkin scheme. We shall denote by x# smooth approximations of \- We 
seek one having the following properties: 

(i) compactly supported Fourier series: = for |k| > K, 

(ii) positivity: \ > 0, 
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FlG. A. 6. Convergence of the fourth- order finite difference scheme, m = 2. Error with respect 
to the exact Dirichlet solution in the interior of the fluid domain (left) and with respect to the the 
penalized solution in the whole domain (right). 



(iii) consistency: \# is close to x m 

(iv) fast decay: is of the order of the numerical round-off error inside f2. 
Condition (i) is necessary for a proper Galerkin discretization, as we shall implement 
further down, while (ii) preserves the positivity of the mollified operator, and (iii) 
keeps the solution close to the one of the original problem. Condition (iv) avoids 
oscillations of x# inside f2, which could create an unacceptable perturbation of the 
solution. 

To enforce (i), we look for as a convolution x# = x*® '■— Jmd x(y ~ x )^Ky)dy, 
where $ is in L 2 (R d ), and the Fourier transform of <&, denoted by $(k), vanishes for 
|k| > K. Condition (ii) implies that $ should be positive, and (iii) means that <E> 
should be as close to 1 as possible. A first idea would be to take $(k) = f if |k| < K 
and otherwise, but that would contradict (ii) and moreover the resulting x # would 
have a very bad localization which enters in conflict with condition (iv). In fact, 
condition (iv) imposes that $ should decay fast at infinity. In [£[, mollifiers having 
optimal localization (in the sense that their variance is minimal in a certain class of 
functions) were constructed as follows. Let Jo be the zeroth order Bessel function of 
the first kind, jo its smallest positive root, and define 

= {'fw fo ; |k| £ 1 (B.D 

10 for |k| > 1. 

By construction, <f> is well localized in the sense that its variance L d |x| 2 0(x)c?x is as 
small as possible, but on the other hand it decays quite slowly at infinity. Indeed, <f> is 
continuous on the circle |k| = I but not C 1 , therefore </>(x) decays at best like |x|" 15 . 
Moreover, </> is not positive. Both issues can be dealt with by taking 

K*\ 2b 



*(x)=C^— I , (B.2) 
where b > I and C is a normalization constant. Such a $ satifies all the conditions 
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FlG. B.l. Comparison of Bessel and Gaussian mollifiers, in physical space (left) and Fourier 
space (right). In physical space, the horizontal scale is expressed in number of grid points. 



that we have imposed, as long as b is chosen properly. In the following, we call <£> a 
Bessel mollifier. 

To choose b we resort to a qualitative judgement, based on Fig. |B.1[ which shows 
the profiles of $ and <£> for b — 1 and 6 = 4. For simplicity K has been scaled to 1. 
We see that b = 4 is a reasonable choice, since <£> decays quickly below the round- 
off error, and the Bessel mollifier with b — 4 roughly extends over 18 grid points. 
For comparison, Gaussian mollifiers are also shown, for which $(k) = cxp (— i<7 2 k 2 ). 
Gaussian mollifiers can be employed, but they are not exactly compactly supported 
in Fourier space, and moreover their localization in physical space is not as good as 
Bessel mollifiers for a given cut-off wavenumber. 
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